Abrupt transition in a sandpile model 
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We present a fixed energy sandpile (FES) model which, by increasing the initial energy, undergoes, 
at the level of individual configurations, a discontinuous transition. The model is obtained by 
modifying the toppling procedure in the BTW Q rules: the energy transfer from a toppling site 
takes place only to neighbouring sites with less energy (negative gradient constraint) and with 
a time ordering (asynchronous). The model is minimal in the sense that removing either of the 
two above mentioned constraints (negative gradient or time ordering) the abrupt transition goes 
over to a continuous transition as in the usual BTW case. Therefore the proposed model offers an 
unique possibility to explore at the microscopic level the basic mechanisms underlying discontinuous 
transitions. 

PACS numbers: 45.70.-n,05.65.+b,05.70.Fh 



INTRODUCTION 



O Self-organized criticality (SOC) has been proposed as a universal mechanism leading to scaling laws in the dynamics 
of driven systems which evolute towards a critical state without the fine tuning of a control parameter pj . Recently it 
has been realized 0, HI 0> IS IE that the approach to criticality is in fact controled by the driving force as well as the 
dissipation present in the dynamics. These observations initiated an extensive search for the influence of the values of 
these control quantities to the critical behaviour of the corresponding SOC models. As a consequence an alternative 
way to study SOC models emerged restricting the considerations in closed systems (without external driving and 
1 — dissipation) obeying the same dynamical rules. In these models the energy density p is exactly conserved (FES). Since 
the fixed-energy sandpiles posses simpler dynamics (without loss or addition) and are translation invariant are easier 
to study than their SOC counterparts. Furthermore the order parameter in this case can be easily identified: below 
a critical density the dynamics lead to an absorbing state characterized by the absense of activity (energy transfer 
t^J- ■ processes) . For densities above the critical p c the system sustains activity. Thus the critical properties of the system 
can be explored defining as order parameter the density of active sites and studying its dependence on p. In addition 
the measured critical exponents in absorbing-state phase transitions can be related to avalanche exponents measured 
' in slowly driven systems 0, Q . 

Our analysis will be devoted exclusively to the sandpile models although many of our results could apply for other 
systems possessing SOC too fOj. Usually the transition from an absor bing (vanishing order parameter) to an active 
state (order parameter different from zero) in FES models is continuous allowing their classification in universality 
classes. On the other hand there is a variety of physical processes which are characterized by a discontinuous transition 
(melting, boiling, earthquake events, etc.). It is therefore natural to ask if such systems could be described in terms 
. of fixed-energy sandpile models. In the present work we show that a suitable modification of the toppling rules in the 
BTW model can lead to a discontinuous transition. Evenmore we determine two conditions, imposed on the BTW 
dynamics, which are both neccessary and sufficient to obtain an abrupt (first order) transition in the corresponding 
FES model: the toppling of an energetically activated site involves energy transfer processes only to less energetic 
neighbours and takes place sequentially in time. This observation opens the possibility to design devices with extreme 
sensitivity on control parameters by applying the analogous constraints. 

Discontinuous transitions in SOC models have already been considered in the context of the so called stick-slip 
dynamics |12| or the breakdown of disordered media |l3j |. However, our approach here is different as it is based 
entirely on slight modifications of the BTW rules. The paper is organized as follows: in Section 2 we present the 
dynamics and describe the critical properties of the proposed model. In addition we compare the obtained results 
with the critical behaviour of conventional FES models. In Section 3 an intepretation (also in terms of microscopic 
dynamics) of the numerical results of Section 2 is given and finally in Section 4 we summarize our findings and give 
a brief outlook. 
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FIG. 1: The function Q a (m) for a single configuration for both the SMBTW FES model using p = 2.0496 as well as the BTW 
FES model (p = 2.1151). 



The FES model which we use is the following: we randomly depose energy on a square lattice in the form of grains. 
Each site can have an arbitrary number of grains. The total energy (and therefore the total number of grains) is fixed 
by the given energy density p. We denote as Zij the number of grains on the site A site is characterized as 

active if Zij exceeds or is equal to a threshold value z c . An active site topples and grains are transfered from this 
site {{i,j)) to the next neighbouring sites provided the corresponding Zi±ij, z^j±x are less than Zij. We use the 
same threshold value z c as in the original BTW model (z c — 4) . We impose however the constraint that the toppling 
procedure is accomplished sequentially: the site (i, j) can topple only if the, suitably defined, previous site has already 
toppled. The ordering of the toppling times of the lattice sites can be chosen randomly with the constraint that the 
entire lattice is covered once. This sequence needs not to be the same in the following sweeps. 

We have explored the dynamics of the system using aLxL square lattice and periodic boundary conditions. As 
usually is the case for the FES models the control parameter is the density of grains p which is a conserved quantity. 
As initial conditions we have used a " microcanonical" ensemble of TV configurations obtained by placing randomly, 
with a uniform distribution, pL 2 grains on the lattice. As an observable to characterize the evolution of the system 
within the ensemble we use the mean density of active sites Q a { i rn) at time m. A stationary state is described through 
the corresponding value of Q — Q a (oo). For most of our numerical simulations we have used N = 1000 and L = 120. 
We observe that for low densities (p < 2.0495) all initial configurations lead to the absorbing state Q — 0. For 
p > 2.0495 (in steps 10 -4 for 100 < L < 316) there is a percentage w of configurations which lead to Q — and the 
remaining configurations lead to Q ~ 0.17. In fact, in the latter case, Q a (m) approaches, in the asymptotic limit (for 
m >> 1), a stationary state described by random fluctuations of amplitude 0.01 around the value 0.17. In Fig. 1 we 
show the function Q a (m) for one such configuration at p = 2.0496. For comparison we display in the same plot the 
corresponding evolution for a typical configuration with density just above the critical one (p = 2.1151) in the usual 
BTW FES model. 

It is clearly seen that there is an energy gap between the stationary states of the system in the case of the SMBTW 
FES model. No such gap can be observed in BTW FES dynamics at this scale. The value of Q just above the critical 
density is two orders of magnitude smaller in the BTW than in the SMBTW case. It is then natural to assign two 
phases to the SMBTW system: phase A is characterized by Q = (absorbing state) while phase B corresponds to 
Q « 0.17. Then the value p c — 2.0495 represents the critical density of the system for the given lattice size and set 
of initial configurations, above which the state with Q ss 0.17 is accessible by the dynamics. In the neighbourhood of 
the critical value the function Q(p) is well fitted by a sigmoidal leading to the estimation of the critical density of the 
SMBTW FES p c , smbtw = 2.0495 ± 0.0072. The set of configurations leading to the nonvanishing value of Q can be 
used to present the phase diagram of the model in the (Q,p) space. As we show in Fig. 2a, at p = p c , smbtw j an 
abrupt jump in Q, possessing the caharacteristics of a first order phase transition, occurs. The plot Q = Q(p) for the 
common BTW model is shown in Fig. 2b. As mentioned in the literature the nature of the phase transition in 



II. THE SMBTW FES SANDPILE MODEL AND ITS CRITICAL BEHAVIOUR 
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FIG. 2: (a) The order parameter Q as a function of the control parameter p of the SMBTW FES model. For p > p c ,smbtw 
we use, in the ensemble averaging, only configurations leading to Q / 0. (b) The corresponding plot for the order parameter 
Q of the FES BTW model. As in (a) we use in the ensemble averaging, for p > p c ,btw, only configurations leading to Q ^ 0. 



this case is not clear due to the devil's staircase form of the function Q{p) [jj, UjJ. Our analysis however supports 
the scenario of a continuous transition in the BTW FES model. This is due to the fact that the spectrum of the 
stationary states P(Q) accessible by the dynamics in the asymptotic limit posesses no energy gap in contrast to the 
SMBTW model. 

In order to support further our conjecture concerning the abrupt behaviour in the SMBTW FES model in contrast 
with a smooth transition of the usual BTW FES model one has to calculate the dependence of the gap of Q on 
the lattice size L. We have performed a calculation of the gap in the SMBTW model for 5 different lattice sizes: 
80, 120, 200, 600, 1000. For lattices with size greater than L — 120 the value of the gap is almost constant: Q « 0.169. 
The results of this calculation are shown with crosses in Fig. 3. A solid line at Q = 0.169 is drawn to guide the 
eye. Additionally we have calculated the asymptotic value Q for the BTW FES model using the same lattice sizes as 
for the SMBTW case. The corresponding results are presented in Fig. 3 with full circles. It is clearly seen that the 
asymptotic value Q in this case tends rapidly to zero with increasing lattice size. Thus, in the thermodynamic limit, 
the gap in the SMBTW remains finite and large, characteristic for a first order transition, while in the BTW case no 
such gap occurs leading to a continuous transition. 

To illustrate this property more transparently we present in Fig. 4 the distribution P(Q) for the two models. To 
allow for a comparison, as the critical density is different in the two cases, we calculate P(Q) at densities p L = p c .i + Sp 
(i = 1,2) using the same value dp for the two models. 



III. INTERPRETATION OF THE NUMERICAL RESULTS 



In order to understand the origin of the sharp transition in the SMBTW model we investigate the corresponding 
dynamics at the microscopic level. We define the single-site states of the system in terms of the possible values of the 
occupancies (number of grains on the site) n s at a given time. To simplify our analysis we include in the state n s = 4 
also the (rare) case when a site is occupied with more that 4 grains. For densities p > p c we observe that the single site 
dynamics are characterized by an ergodic behaviour: each one-site trajectory visits irregularly all the accesible states 
in phase space. After a characteristic time scale t w 2000 (algorithmic time in units of lattice sweeps) an invariant 
density is established. Having achieved this stationarity each state is visited with almost equal probability by the 
dynamically evolving site. A very smooth maximum occurs for n s — 1 and n s = 3. As the number of possible states 
is 5 (n s — 0, 1, 2, 3, 4) we expect, for a uniform invariant density and assuming that the ergodic hypothesis applies, to 
have 1/5 probability to be on the active state (n s = 4), a value leading to Q — 0.2 very close to the observed value 
Q * 0.17. The deviation is due to the fact that the invariant density is not exactlty uniform. 

It is worth exploring the global dynamics of the system as well. Therefore we investigate the evolution of a typical 



4 



o 



10" 1 



10"' 



10"' 



10"' 



• SMBTW 
A BTW 



200 400 600 800 1000 

L 



FIG. 3: The gap in Q as a function of the lattice size L for the SMBTW (full circles) as well as the BTW (full triangles) FES 
model (semilog plot). The solid line at Q = 0.169 is used to guide the eye. 




FIG. 4: The distribution P(Q) for the SMBTW and BTW models at densities pi = 2.060 and p 2 = 2.1255 respectively 
(5p = 0.0105 for both models). 
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FIG. 5: The evolution of a typical configuration using the SMBTW rules for p = 2.060. We present the initial state (a) as 
well as the resulting state for t = 5000 (b). As explained in the text only a part of the lattice is displayed. 

configuration of the entire lattice in the critical region p > p c . In Fig. 5 we present the evolution of such a configuration 
for p — 2.060. We show the contourplots for the initial configuration as well as the envolving state at algorithmic time 
t = 5000. We recall here that stationarity is achieved already at t = 2000. To simplify our representation without 
loosing on physical insight we adopt the coarsegrained description of the phase space of the system used in 0, ITt^ 
distinguishing between stable (n s < 3, gray), critical (n s — 3, white) and active (n s > 3, black) sites. For a better 
visualization of the details of the dynamics it is convenient to zoom into a part of the lattice which we choose here to 
be the set of sites with i,j G [40,60]. Obviously, as displayed in Fig. 5a, initially the active sites form clusters 
which are embedded in domains consisting of critical sites. 

As the dynamics evolve the active clusters split and their gradual destruction initiates in favour of the formation 
of larger domains of critical sites. Before the active clusters disappear the domains of critical sites approaching 
each other collide forming larger domains of critical sites containing smaller active clusters. When the stationary 
state is achieved the critical sites form large connected domains, similar to the above-threshold percolation clusters, 
which span the entire lattice. A dynamical equilibrium is established characterized by an irregular deformation of 
the critical domains as well as a chaotic motion, icorporating the spliting and recombining through collisions, of the 
active clusters within these domains. Thus the critical domains construct a protecting network of communicating 
channels (see Fig. 5b) for the irregular evolution of the active clusters. It must be noted that critical clusters can be 
formed also in the conventional sandpile models. However the maximal size of critical clusters in the case of the FES 
SMBTW model is much larger than the corresponding size of the conventional models. 

The formation of this network is a process depending very sensitively on the initial conditions. Consider for example 
a typical configuration leading to Q — for p — 2.0495. The total number of grains, for a L = 120 lattice, in such 
a configuration is 25912. Inserting one additional grain in the lattice can lead (depending on where we put it) to 
Q w 0.17. Thus, density fluctuations of the order of 10~ 5 may lead to a tremendous change of the order parameters' 
value. We attribute this behaviour to the first order character of the transition in our model. Actually imposing a 
time ordering in the toppling rules of the system introduces an internal timescale of the order of the lattice size L. 
Changes in the environment of each site during time intervals of this order are felt by the corresponding site due to 
the dependence of the local energy current on the configuration of the environment expressed through the negative 
gradient constraint. This feedback mechanism creates an unstable environment which leads to a strong sensitivity to 
initial conditions. The relative timescale for which this instability influences the dynamics is of the order of 4- and 
therefore environmental influence becomes continuous in time in the limit of an infinite system. We have investigated 
the transition by relaxing each time one of the two constraints we have used in the SMBTW model, i.e. the negative 
gradient and the time ordering. In both cases the transition turns out to be continuous and the corresponding phase 
diagram is very close to the diagram shown in Fig. 2b for the BTW model. Keeping only the time ordering constraint 
(SBTW model) the internal timescale has no consequence on the evolution of the system as the corresponding energy 
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currents do not depend on the environment. Keeping only the negative gradient constraint (PMBTW model) there is 
no internal time and the changes in the environment occur in timescales of the order of L 2 which become less and less 
important as the system size increases. On the other hand the choice of negative gradient (instead of, for example, a 
positive gradient) constraint is neccessary in order to reach a stationary absorbing state. 

IV. CONCLUDING REMARKS 

Let us now summarize briefly our results. In the present work we have introduced a sandpile model resulting from 
the BTW rules through the addition of two constraints: energy is transfered from an active site only to less energetic 
neighbouring sites and the toppling takes place in an asynchronous manner. The later means that the instants when 
a lattice site can topple are ordered in time. We have investigated the case without external driving. The system 
undergoes an abrupt (first order) transition. Our model resembles the dynamics of activated random walkers 01 
with the additional property of walking only to less occupied sites. This introduces a feedback mechanism influencing 
locally the energy flow on the lattice. The proposed model is minimal in the sense that both constraints are neccessary 
in order to achieve the discontinuous change of the order parameter. Considering the one-site as well as the global 
dynamics of the system we were able to understand the qualitative as well as some quantitative features of the basic 
mechanism underlying the observed transition. Our model is capable to describe physical systems characterized 
by asynchronous spontaneous energy transfer. In particular the proposed model could give additional insight in 
prefracture processes during earthquake (EQ) events. In this case the imposed two constraints are fulfiled: (i) the 
stress within the focal area is transfered from regions of higher tension to regions of lower tension (negative gradient) 
while (ii) the microcrack transmission takes place sequentially. When the fracture of the heterogeneous environment 
is consumated the remaining asperities suffer an intensive stress from their environment. A tiny fluctuation of the 
surrounding stress field is decisive for the final fracture of asperities and therefore the occurance of an EQ event or not. 
In the former case the corresponding transition is abrupt (first order) |l8t Il9| . Additionally the extreme sensitivity 
of the proposed model on density fluctuations of relative magnitude of the order of 10 -5 (or even less for lattices 
with L > 120) suggests the possibility to design a high efficiency sensor by appropriate realization of the proposed 
dynamical rules. It remains a challenging task to determine in what extent the observed behaviour possesses universal 
features. 
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